#######################################################
####                                               ####
####  Input: MP-by-debate data for analysis        ####
####  Output: Supplementary materials analysis     ####
####                                               ####
#######################################################

stdize <- function(x) (x - mean(x, na.rm = T))/sd(x, na.rm = T)

# Load libraries

library(data.table) # v1.14.8
library(plyr) # v1.8.9
library(ggplot2) # v3.4.3
library(sandwich) # v3.0-2
library(quanteda) # v3.3.1
library(dplyr) # v1.1.3
library(sjPlot) # v2.8.15
library(sjmisc) # v2.8.9
library(broom) # v1.0.5
library(plm) # v2.6-3
library(estimatr) #v1.0.0

# Load data

load("working/speech_scores.Rdata")

# Regression --------------------------------------------------------------------------------------------------------------

women_1 <-
  lm(women_dic_std ~ gender * years_in_parliament, data = speech_scores)
women_2 <-
  lm(
    women_dic_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
women_3 <-
  lm(
    women_dic_std ~ gender * years_in_parliament + party + age_years + women_minister_gov + gov_mp +
      women_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
women_4 <- plm(
  women_dic_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
women_5 <-
  plm(
    women_dic_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

women_4_se <- sqrt(diag(vcovCR(
  women_4,
  cluster = "individual",
  type = "CR0"
)))

women_5_se <- sqrt(diag(vcovCR(
  women_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(women_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <- sqrt(diag(vcovCL(women_2, ~ person_id)))
my_clustered_errors[[3]] <- sqrt(diag(vcovCL(women_3, ~ person_id)))
my_clustered_errors[[4]] <- women_4_se
my_clustered_errors[[5]] <- women_5_se

save(women_1,
     women_2,
     women_3,
     women_4,
     women_5,
     my_clustered_errors,
     file = "analysis/models/women.Rdata")

health_1 <-
  lm(health_std ~ gender * years_in_parliament, data = speech_scores)
health_2 <-
  lm(
    health_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
health_3 <-
  lm(
    health_std ~ gender * years_in_parliament + party + age_years + health_minister_gov + gov_mp +
      health_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
health_4 <- plm(
  health_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
health_5 <-
  plm(
    health_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

health_4_se <- sqrt(diag(vcovCR(
  health_4,
  cluster = "individual",
  type = "CR0"
)))

health_5_se <- sqrt(diag(vcovCR(
  health_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(health_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <- sqrt(diag(vcovCL(health_2, ~ person_id)))
my_clustered_errors[[3]] <- sqrt(diag(vcovCL(health_3, ~ person_id)))
my_clustered_errors[[4]] <- health_4_se
my_clustered_errors[[5]] <- health_5_se

save(health_1,
     health_2,
     health_3,
     health_4,
     health_5,
     my_clustered_errors,
     file = "analysis/models/health.Rdata")

education_1 <-
  lm(education_std ~ gender * years_in_parliament, data = speech_scores)
education_2 <-
  lm(
    education_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
education_3 <-
  lm(
    education_std ~ gender * years_in_parliament + party + age_years + education_minister_gov + gov_mp +
      education_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
education_4 <- plm(
  education_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)

education_5 <-
  plm(
    education_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

education_4_se <- sqrt(diag(vcovCR(
  education_4,
  cluster = "individual",
  type = "CR0"
)))

education_5_se <- sqrt(diag(vcovCR(
  education_5,
  cluster = "individual",
  type = "CR0"
)))
my_clustered_errors <-
  starprep(education_1,
           clusters = speech_scores$person_id,
           se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(education_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(education_3, ~ person_id)))
my_clustered_errors[[4]] <- education_4_se
my_clustered_errors[[5]] <- education_5_se

save(
  education_1,
  education_2,
  education_3,
  education_4,
  education_5,
  my_clustered_errors,
  file = "analysis/models/education.Rdata"
)

social_1 <-
  lm(social_std ~ gender * years_in_parliament, data = speech_scores)
social_2 <-
  lm(
    social_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
social_3 <-
  lm(
    social_std ~ gender * years_in_parliament + party + age_years + welfare_minister_gov + gov_mp +
      welfare_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
social_4 <- plm(
  social_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
social_5 <-
  plm(
    social_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

social_4_se <- sqrt(diag(vcovCR(
  social_4,
  cluster = "individual",
  type = "CR0"
)))
social_5_se <- sqrt(diag(vcovCR(
  social_5,
  cluster = "individual",
  type = "CR0"
)))
my_clustered_errors <-
  starprep(social_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <- sqrt(diag(vcovCL(social_2, ~ person_id)))
my_clustered_errors[[3]] <- sqrt(diag(vcovCL(social_3, ~ person_id)))
my_clustered_errors[[4]] <- social_4_se
my_clustered_errors[[5]] <- social_5_se

save(social_1,
     social_2,
     social_3,
     social_4,
     social_5,
     my_clustered_errors,
     file = "analysis/models/social.Rdata")

environment_1 <-
  lm(environment_std ~ gender * years_in_parliament, data = speech_scores)
environment_2 <-
  lm(
    environment_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
environment_3 <-
  lm(
    environment_std ~ gender * years_in_parliament + party + age_years + environment_minister_gov + gov_mp +
      environment_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
environment_4 <- plm(
  environment_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
environment_5 <-
  plm(
    environment_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

environment_4_se <- sqrt(diag(vcovCR(
  environment_4,
  cluster = "individual",
  type = "CR0"
)))

environment_5_se <- sqrt(diag(vcovCR(
  environment_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(environment_1,
           clusters = speech_scores$person_id,
           se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(environment_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(environment_3, ~ person_id)))
my_clustered_errors[[4]] <- environment_4_se
my_clustered_errors[[5]] <- environment_5_se

save(
  environment_1,
  environment_2,
  environment_3,
  environment_4,
  environment_5,
  my_clustered_errors,
  file = "analysis/models/environment.Rdata"
)

defence_1 <-
  lm(defence_std ~ gender * years_in_parliament, data = speech_scores)
defence_2 <-
  lm(
    defence_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
defence_3 <-
  lm(
    defence_std ~ gender * years_in_parliament + party + age_years + defence_minister_gov + gov_mp +
      defence_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
defence_4 <- plm(
  defence_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
defence_5 <-
  plm(
    defence_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )


defence_4_se <- sqrt(diag(vcovCR(
  defence_4,
  cluster = "individual",
  type = "CR0"
)))

defence_5_se <- sqrt(diag(vcovCR(
  defence_5,
  cluster = "individual",
  type = "CR0"
)))
my_clustered_errors <-
  starprep(defence_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(defence_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(defence_3, ~ person_id)))
my_clustered_errors[[4]] <- defence_4_se
my_clustered_errors[[5]] <- defence_5_se

save(defence_1,
     defence_2,
     defence_3,
     defence_4,
     defence_5,
     my_clustered_errors,
     file = "analysis/models/defence.Rdata")

economy_1 <-
  lm(economy_std ~ gender * years_in_parliament, data = speech_scores)
economy_2 <-
  lm(
    economy_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
economy_3 <-
  lm(
    economy_std ~ gender * years_in_parliament + party + age_years + economy_minister_gov + gov_mp +
      economy_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
economy_4 <- plm(
  economy_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
economy_5 <-
  plm(
    economy_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

economy_4_se <- sqrt(diag(vcovCR(
  economy_4,
  cluster = "individual",
  type = "CR0"
)))

economy_5_se <- sqrt(diag(vcovCR(
  economy_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(economy_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(economy_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(economy_3, ~ person_id)))
my_clustered_errors[[4]] <- economy_4_se
my_clustered_errors[[5]] <- economy_5_se

save(economy_1,
     economy_2,
     economy_3,
     economy_4,
     economy_5,
     my_clustered_errors,
     file = "analysis/models/economy.Rdata")

agriculture_1 <-
  lm(agriculture_std ~ gender * years_in_parliament, data = speech_scores)
agriculture_2 <-
  lm(
    agriculture_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
agriculture_3 <-
  lm(
    agriculture_std ~ gender * years_in_parliament + party + age_years + agriculture_minister_gov + gov_mp +
      agriculture_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
agriculture_4 <- plm(
  agriculture_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
agriculture_5 <-
  plm(
    agriculture_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )

agriculture_4_se <- sqrt(diag(vcovCR(
  agriculture_4,
  cluster = "individual",
  type = "CR0"
)))
agriculture_5_se <- sqrt(diag(vcovCR(
  agriculture_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(agriculture_1,
           clusters = speech_scores$person_id,
           se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(agriculture_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(agriculture_3, ~ person_id)))
my_clustered_errors[[4]] <- agriculture_4_se
my_clustered_errors[[5]] <- agriculture_5_se

save(
  agriculture_1,
  agriculture_2,
  agriculture_3,
  agriculture_4,
  agriculture_5,
  my_clustered_errors,
  file = "analysis/models/agriculture.Rdata"
)

trade_1 <-
  lm(trade_std ~ gender * years_in_parliament, data = speech_scores)
trade_2 <-
  lm(
    trade_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
trade_3 <-
  lm(
    trade_std ~ gender * years_in_parliament + party + age_years + trade_minister_gov + gov_mp +
      trade_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
trade_4 <- plm(
  trade_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
trade_5 <-
  plm(
    trade_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )
trade_4_se <- sqrt(diag(vcovCR(
  trade_4,
  cluster = "individual",
  type = "CR0"
)))
trade_5_se <- sqrt(diag(vcovCR(
  trade_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(trade_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <- sqrt(diag(vcovCL(trade_2, ~ person_id)))
my_clustered_errors[[3]] <- sqrt(diag(vcovCL(trade_3, ~ person_id)))
my_clustered_errors[[4]] <- trade_4_se
my_clustered_errors[[5]] <- trade_5_se

save(trade_1,
     trade_2,
     trade_3,
     trade_4,
     trade_5,
     my_clustered_errors,
     file = "analysis/models/trade.Rdata")

crime_1 <-
  lm(crime_std ~ gender * years_in_parliament, data = speech_scores)
crime_2 <-
  lm(
    crime_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
crime_3 <-
  lm(
    crime_std ~ gender * years_in_parliament + party + age_years + crime_minister_gov + gov_mp +
      crime_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
crime_4 <- plm(
  crime_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
crime_5 <-
  plm(
    crime_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )
crime_4_se <- sqrt(diag(vcovCR(
  crime_4,
  cluster = "individual",
  type = "CR0"
)))
crime_5_se <- sqrt(diag(vcovCR(
  crime_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(crime_1, clusters = speech_scores$person_id, se_type = "stata")
my_clustered_errors[[2]] <- sqrt(diag(vcovCL(crime_2, ~ person_id)))
my_clustered_errors[[3]] <- sqrt(diag(vcovCL(crime_3, ~ person_id)))
my_clustered_errors[[4]] <- crime_4_se
my_clustered_errors[[5]] <- crime_5_se

save(crime_1,
     crime_2,
     crime_3,
     crime_4,
     crime_5,
     my_clustered_errors,
     file = "analysis/models/crime.Rdata")

transport_1 <-
  lm(transport_std ~ gender * years_in_parliament, data = speech_scores)
transport_2 <-
  lm(
    transport_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
transport_3 <-
  lm(
    transport_std ~ gender * years_in_parliament + party + age_years + transport_minister_gov + gov_mp +
      transport_minister_opp + marginality + occupation_class,
    data = speech_scores
  )
transport_4 <- plm(
  transport_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
transport_5 <-
  plm(
    transport_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )
transport_4_se <- sqrt(diag(vcovCR(
  transport_4,
  cluster = "individual",
  type = "CR0"
)))
transport_5_se <- sqrt(diag(vcovCR(
  transport_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(transport_1,
           clusters = speech_scores$person_id,
           se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(transport_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(transport_3, ~ person_id)))
my_clustered_errors[[4]] <- transport_4_se
my_clustered_errors[[5]] <- transport_5_se

save(
  transport_1,
  transport_2,
  transport_3,
  transport_4,
  transport_5,
  my_clustered_errors,
  file = "analysis/models/transport.Rdata"
)

children_1 <-
  lm(children_std ~ gender * years_in_parliament, data = speech_scores)
children_2 <-
  lm(
    children_std ~ gender * years_in_parliament + party + age_years + attends_cabinet +
      attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
      is_committee_chair + marginality + occupation_class,
    data = speech_scores
  )
children_3 <-
  lm(
    children_std ~ gender * years_in_parliament + family_minister_opp + marginality + gov_mp
    + occupation_class,
    data = speech_scores
  )
children_4 <- plm(
  children_std ~ gender * years_in_parliament,
  index = "person_id",
  model = "within",
  effect = "individual",
  data = speech_scores
)
children_5 <-
  plm(
    children_std ~ gender * years_in_parliament + factor(parliamentary_term),
    index = "person_id",
    model = "within",
    effect = "individual",
    data = speech_scores
  )
children_4_se <- sqrt(diag(vcovCR(
  children_4,
  cluster = "individual",
  type = "CR0"
)))
children_5_se <- sqrt(diag(vcovCR(
  children_5,
  cluster = "individual",
  type = "CR0"
)))

my_clustered_errors <-
  starprep(children_1,
           clusters = speech_scores$person_id,
           se_type = "stata")
my_clustered_errors[[2]] <-
  sqrt(diag(vcovCL(children_2, ~ person_id)))
my_clustered_errors[[3]] <-
  sqrt(diag(vcovCL(children_3, ~ person_id)))
my_clustered_errors[[4]] <- children_4_se
my_clustered_errors[[5]] <- children_5_se

save(
  children_1,
  children_2,
  children_3,
  children_4,
  children_5,
  my_clustered_errors,
  file = "analysis/models/children.Rdata"
)

## Covariate effects --------------------------------------------------------------------------------------------------------------

speech_scores$age_std <- stdize(speech_scores$age_years)

children <-
  lm_robust(
    children_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

health <-
  lm_robust(
    health_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

education <-
  lm_robust(
    education_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

women <-
  lm_robust(
    women_dic_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

social <-
  lm_robust(
    social_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

crime <-
  lm_robust(
    crime_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

environment <-
  lm_robust(
    environment_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

transport <-
  lm_robust(
    transport_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

economy <-
  lm_robust(
    economy_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

trade <-
  lm_robust(
    trade_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

agriculture <-
  lm_robust(
    agriculture_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

defence <-
  lm_robust(
    defence_std ~ gender + party + age_std + attends_cabinet + attends_shadow_cabinet + is_gov_minister +
      is_opp_minister + is_committee_chair + marginality + years_in_parliament + gov_mp + occupation_class,
    data = speech_scores,
    se_type = "stata",
    clusters = person_id
  )

covariates_for_plot <- bind_rows(
  tidy(children, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Children & Family"),
  tidy(health, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Health"),
  tidy(education, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Education"),
  tidy(women, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Women"),
  tidy(social, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Welfare"),
  tidy(social, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Crime & Policing"),
  tidy(environment, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Environment"),
  tidy(transport, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Transport"),
  tidy(economy, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Finance & Economy"),
  tidy(trade, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Trade"),
  tidy(agriculture, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Agriculture"),
  tidy(defence, conf.int = T)[c(2:19), ] %>% mutate(outcome = "Defence"),
)

covariates_for_plot$sig <-
  ifelse(
    abs(covariates_for_plot$estimate / covariates_for_plot$std.error) >= 1.96,
    TRUE,
    FALSE
  )

pdf("analysis/plots/figure_S7.pdf", 11.5, 7)
ggplot(covariates_for_plot, aes(y = term, x = estimate)) +
  geom_point(color = ifelse(covariates_for_plot$sig, "black", alpha("black", .2))) +
  ylab("") + xlab("") +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high),
    color = ifelse(covariates_for_plot$sig, "black", alpha("black", .2))
  ) +
  facet_wrap( ~ outcome, scales = "free") +
  scale_y_discrete(
    labels = c(
      "Age",
      "Cabinet",
      "Shadow Cabinet",
      "Woman MP",
      "Government MP",
      "Com. Chair",
      "Gov. Minister",
      "Opp. Minister",
      "Margin of victory",
      "Occ.: Professional",
      "Occ.: Political",
      "Occ.: Other",
      "Occ.: Manual",
      "Party: Green",
      "Party: Labour",
      "Party: Lib Dem",
      "Party: Other",
      "Years in Parliament"
    )
  )
dev.off()

# Talking about women within debates on all policy areas

children <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                 attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                 is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$children_debate==TRUE])
health <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
               attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
               is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$health_debate==TRUE])
education <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                  attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                  is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$education_debate==TRUE])
welfare <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$welfare_debate==TRUE])
crime <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
              attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
              is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$crime_debate==TRUE])
environment <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                    attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                    is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$environment_debate==TRUE])
transport <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                  attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                  is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$transport_debate==TRUE])
economy <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$economy_debate==TRUE])
trade <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
              attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
              is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$trade_debate==TRUE])
agriculture <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                    attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                    is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$agriculture_debate==TRUE])
defence <- lm(women_dic_std ~ gender*years_in_parliament + party + attends_cabinet +
                attends_shadow_cabinet + is_gov_minister + is_opp_minister + gov_mp +
                is_committee_chair + marginality + occupation_class, data = speech_scores[speech_scores$defence_debate==TRUE])

save(children, health, education, welfare, crime, environment, transport, economy, trade, agriculture, defence, file = "analysis/models/debates_seniority.Rdata")


